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Abstract 

We develop a physiologically-based lattice model for the transport and metabolism 
of drugs in the functional unit of the liver, called the lobule. In contrast to earlier 
studies, we have emphasized the dominant role of convection in well-vascularized 
tissue with a given structure. Estimates of convective, diffusive and reaction 
contributions are given. We have compared drug concentration levels observed 
exiting the lobule with their predicted detailed distribution inside the lobule, 
assuming that most often the former is accessible information while the latter is not. 
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Background 

The liver is the major site of biotransformation of endogenous and xenobiotic 
substances including drugs in the body. Its main role is to prevent accumulation of a 
wide range of chemical compounds in the blood by converting them into a form 
suitable for elimination. Such vital processes, however, can potentially damage liver 
tissue and hence its functionality. Examinations of hepatic clearance shows that sub- 
stance extraction not only can be limited by damaged hepatocytes (liver cells) but also 
by the intrinsic (enzymatic) ability to eliminate the drug, by the resistance to drug 
transport from blood to eliminating tissue cells, and by the hepatic blood flow itself. 
Indeed, perturbations in the hepatic flow patterns e.g. through disease or aging can 
significantly alter the systemic clearance of these substances. As a result, a quantitative 
understanding between the liver performance and its structural integrity would be a 
great utility in the toxicology of newly developed drugs. 

To quantitatively assess these interacting processes, numerous models of hepatic 
clearance have been developed both in vitro and in silico, see an extensive review in 
[1]. Although in vitro liver models can be considered as the best alternative to analyze 
actual organ performance, their major shortcoming is to assume the liver as a homo- 
geneous environment. It is known that metabolic functions, such as xenobiotic meta- 
bolism, amino acid conversion, cholesterol synthesis, oxidation, etc. are all zone 
specific. As a result, in vitro models are not necessarily well suited for studying and 
capturing the physiology of liver with hepatocyte morphology, property and functio- 
nality varying from point to point throughout of the organ. Because of this, there has 
been surge of interest among investigators in developing computational (in silico) 
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models of the liver in order to study inhomogeneous structures in the presence of 
xenobiotic metabolisms and their effects on the hepatic clearance. A major challenge 
remains to integrate both physiological aspects and various metabolic processes in a 
single model. 

Generally, in all computational models, nutrient and metabolite containing blood 
enters through a segment of interest (the lobule) via portal veins (and hepatic arteries) 
and flows through small channels (sinusoids) that circumvent liver cells (hepatocytes), 
Figure la. These (possibly toxic) chemicals transport from sinusoids to tissue, are 
converted to other forms by hepatocyte metabolic reactions, transport back to the 
sinusoids, and are finally excreted from the lobule via the central (hepatic) vein. 

Multiple computational models have been proposed: single and multi-compartment 
models, distribution-based models, agent-based models, interconnected parallel tubes, 
etc. See [1] for further details. Two simple and commonly employed models are (a) the 
well-stirred compartment model, and (b) the parallel tube model. The first considers 
the liver as a homogeneous compartment [2] and the concentration of drug in the liver 
to be in equilibrium with venous (emerging) blood drug concentration. The second, 




(c) (d) 

Figure 1 Flow network structure, (a) Schematic diagram of a cross section of hepatic parenchyma 

consisting hexagonal lobules, portal and hepatic veins. The lobule contains liver cells (sinusoids and 

hepatocytes). The segment represents a typical area studied in this paper, (b) Homogeneous lattice 

(segmented area) with high porosity bands (in red) representing sinusoids and lower porosity regions (in 

blue) representing tissue containing hepatocytes. (c) 2D hexagonal lobule lattice with sinusoids generated 

via a diffusion limited aggregation algorithm, (d) 3D hexagonal lobule lattice. 
\ . J 
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developed by Bass and coworkers [3,4], regards the liver as a series of parallel tubes 
and the concentration of drug declines along the length of these tubes, due to enzymes 
distributed evenly along the length of these tubes. Pang and Rowland [5,6] present an 
in-depth comparison of the strengths and limitations of each approach. Subsequent 
developments [7-11] indicate that these two models represent limiting cases, and 
attempts to generalize each, e.g. by adding a series of well stirred compartments or by 
considering a distribution of tubes or adding directly dispersive mixing, have tended to 
make the models more consistent with each other. Further enhancements of these 
models include consideration of heterogeneity by adding stochastic processes, 
Michaels-Menten kinetics to describe intrinsic elimination, uniform and non-uniform 
metabolic processes, and enzyme zonation along the sinusoid by stacking a number of 
compartments, each with different concentration level and metabolic activities [1]. We 
emphasize, however, all of the above models are essentially ID or pseudo ID 
approaches and hence ignore the true three dimensional structure of the liver. 

In series of papers, we intend to introduce a physiologically-based lattice model of 
the functional unit of the liver that integrates all structural and physiological aspects 
discussed above. Our model takes into account parameters such as the distribution 
volume, permeability, and porosity of the liver vasculature and cells. Furthermore, it 
includes both flow-limited and diffusion-limited exchange of drug molecules into the 
extravascular space from the sinusoids; sequestration of the drug molecules within liver 
cells with enzymatic transformation; and exchange of the metabolized drug molecules 
back from liver cells to the vasculature. Estimates and consequences of the competing 
flow processes are given. Furthermore, the enzymatic transformation of the drug can 
be either simple or saturable. The model allows us to include the effects of the 
intrahepatic mixing process on the enzymatic transformation of drug molecules at the 
cellular level. 

In this paper, we focus on a simple regular lattice to define a base-case model for our 
study. Here we explore the dynamics of competing convective, diffusive, and reactive 
processes acting on an injected drug. Multiple sensitivity simulations are performed 
and their consequences on drug concentration levels found exiting the lobule as well as 
their detailed spatial distribution within the lobule are discussed. Future extensions to 
additional structural and/or physiological inhomogeneities such as non-regular and 
statistical lattices, enzymatic zonation, and extensions to 2D / 3D hexagonal lattices will 
be developed, see Figure lc and Id. 

The liver and drug kinetics 

Liver architecture 

At the macroscopic level, the liver consists of three vascular trees, two supply trees that 
originate from the portal vein and hepatic artery, and one collecting tree that drains 
into the portal artery [12]. The vessels bifurcate down to the terminal arterioles and ve- 
nules, which are organized into portal tracts along with a terminal bile duct. Liver cells, 
called hepatocytes, radiate outward from the terminal vessels. These plates of hepato- 
cytes are interspersed by sinusoids, which play the role of the capillary in the liver, and 
the spaces of Disse, which are the extravascular space of the liver [13]. Finally, the 
blood is collected and removed by the hepatic venules. 



Rezania et al. Theoretical Biology and Medical Modelling 2013, 10:52 
http://www.tbiomed.eom/content/1 0/1 /52 



Page 4 of 1 9 



Functional unit 

The functional unit of an organ is the smallest structural unit that can independently 
serve all of the organs functions [14]. Because of its complexity, there is continued 
debate about what the functional unit of the liver should be. The classic lobule is a 
hexagonal cylinder, centered around a hepatic venule and with portal tracts situated at 
the corners. The portal lobule has a similar shape but is centered about a portal tract 
with the hepatic venules at the periphery [15]. The acinus is another proposed unit and 
is based on the pattern formed by the cords of hepatocytes between two central 
venules. Matsumoto and Kawakami [16] suggested that the classic lobule can be 
divided into primary lobules, which are cone-shaped and each fed by one portal tract 
and drained by one hepatic venule. Teutsch and colleagues [17,18] performed a mor- 
phological study of rat and human liver lobules, and their results support the idea of a 
secondary unit made up of primary units in what they term as a modular architecture. 
They conclude, however, that the primary unit is more polyhedral in shape than 
conical. Other experiments done by Ruijter et al. [19] suggest that the primary unit is 
needle-shaped and that there are equal amounts of portal and central vein associated 
with one unit. For this study, the primary unit is taken to be one-fourth of the classical 
lobule. The relevant anatomical values are listed in Table 1. 

Elimination kinetics 

Here drug uptake and elimination (i.e. conversion to metabolized product) is viewed as 
a single-step saturable process following Michaelis-Menten kinetics [20], such that 



Here C(t) is the local drug concentration and dC(t)/dt is the drug metabolization 
rate. (In what follows, C(t) is expressed as px t with p the fluid density and x t mole 
fractions of /-th species). Note we are explicitly modeling the drug transport to an indi- 
vidual hepatocyte surface via our lobule lattice model and assume an effective one- step 
reaction transformation beyond that point. We recognize that drug incorporation and 
elimination is still a multi-step process even once the drug reaches the cell surface 
however. It is hoped that these approximations ignore processes occurring on a shorter 
time scale than the experimental resolution. Nonetheless, we explore possible compli- 
cations via simple sensitivities to the choice of reaction time constant. 

Table 1 Anatomical parameter values for the liver [12,13] 

Parameter Value 

Hepatocyte diameter 12 -24 urn 

Diameter of liver cell sheets 25 urn 

Lobule diameter 1 - 2.5 mm 

Mean sinusoid diameter 7.3 urn 

Vascular tissue component 28 - 30 % 

Specific gravity of liver 1 .05 

Liver volume 1 071 ± 228 cm 3 



dC(t) 



max 

C(t)/(K m + C(t)) 



(1) 



dt 
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Model and method 

The primary unit of the liver was approximated by a symmetry element of a 51 x 51 
square lattice such that four units make up one lobule. The architecture of the lattice 
consists of hepatocyte grid cells interlaced by a network of narrower sinusoidal grid 
cells (Figure la). This "intermediate realism" model represents a practical compromise 
with simulation run times of seconds compared with precise, high resolution, high 
realism models which we will explore in future work. These latter models can be 
expected to have equivalent run times of hours and before proceeding to such simula- 
tions, it was deemed prudent to explore a wide range of process variables including 
both lattice effects as well as metabolic effects with our simpler representation in order 
to get an overall assessment of liver lobule performance and drug metabolism in 
healthy and diseased livers. Additionally, realistic, high resolution models by their very 
nature can contain conflicting mechanisms whereas simpler models allow their sepa- 
ration for analysis. 

In our model, the vascular supply tree (hereafter termed "injector") is located at 
top left-corner (grid cell 1,1) while the vascular collecting tree (hereafter termed 
"producer") is located at bottom right-corner (grid cell 51,51) of this model. A 
complete lobule would be formed by placing three more equivalent lattices around 
the central injector but because of symmetry for this ideal lattice model, all 4 sec- 
tions would produce equivalent profiles and hence be computationally redundant. 
The diameter of the sinusoid grid cells was taken to be 0.0006 cm, and the dia- 
meter of the hepatocyte grid cells was taken to be 0.0024 cm. (Again because of 
symmetry, the bounding sinusoids are taken as % size, 0.0003 cm, as the remaining 
portion of these sinusoids are included in the adjacent sections). The length of the 
lattice is thus 0.0750 cm per side. Doubling this value gives a lobule diameter of 
0.150 cm, which is consistent with values listed in Table 1. 

Convective molar flux is modelled according to Darcys Law [21]: 



where f ik is the z-th component of fluid flux in /^-direction, p and ft are fluid molar 
density and viscosity, v k , K k and V/^> are the Darcy velocity, permeability, and pressure 
gradient in direction /<, respectively. The blood viscosity ft is taken to be 3.5 mPa-s 
(3.5 centipoise). Blood molar density is assumed that of water, p at 55.4 mmoles/cm 3 . It 
is emphasized that in this paper, following Darcys Law and the conventions of flow in 
porous media, permeability K is defined as a measure of the transmissibility of a grid 
cell to the flow of a fluid, and is expressed in units of area (e.g. cm 2 ). 

Each sinusoid grid cell represents a tubular vessel of diameter 2a. Taking the ratio of 
the volume of the vessel to the volume of the grid cell yields a porosity of 



J%. = pXiV k = pxi — V k p 



ft 



(2) 



0si: 



71 a 2 (2a) 71 



= - = 0.7854 



(3) 



sin 



{2af 4 
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For a cylindrical tube, the permeability can be calculated from the tube radius [14] 
K sin =^ = 1.126 /mi 2 (4) 

o 

Here and in what follows, it is noted that porosity is dimensionless while permeability 
has units of a characteristic length squared. 

Each parenchymal grid cell represents a cellular (hepatocyte) component and an 
extracellular (space of Disse) component. A ratio of 0.75 to 0.25 was chosen for their 
respective contributions to the volume, and the porosity of the parenchymal sites was 
therefore 0.25. 

The corresponding permeability of the tissue grid cells is estimated from a Carmen- 
Kozeny model [21] of flow around a spherical object (the hepatocyte) of radius 12 |im. 
This states the permeability is proportional to the objects diameter D p and the porosity 
0 tis as follows: 

K tis = . (5) 

180(l-fc) 2 

An ideal result where D p /L = 1 or R/L = 0.5 with L as grid size and R as particle radius 
gives 

(1-0 J = = 0.5236 or <p tis = 0.4764 (6) 

Then assuming as we do, L = D p = 24 (im, the Carmen-Kozeny formula yields 
K tis = 1.262 ftm 2 

The Carman-Kozeny expression basically states physically that the order of magni- 
tude of the permeability scales with the particle size squared. 

The above analysis is an ideal result as the ECM around the cell particle will further 
reduce porosity. Thus we could study a range of tissue effective porosities, leading to a 
range of effective tissue permeabilities. Table 2 summarizes a range of possible values. 
In this paper, we will utilize the "base case" value for tissue porosity and permeability, 
while the effects of more extreme choices will be examined in a second paper. 

The convective driving force originates from an input site corresponding to a ter- 
minal portal venule at one corner of the lattice and an output site corresponding to a 
terminal hepatic venule at the opposite corner. For simplicity, the hepatic artery blood 
supply, which is lower in volume and pulsatile in nature, is omitted for the current sim- 
ulations. The pressure value at the inlet and the outlet are taken to be P in - 103 kPa 
and P out = 101.8 kPa, respectively. After subtracting the atmospheric pressure, these 
values are consistent with experimental values quoted by Rappaport [22], who found 
that the terminal portal venule pressure was in the range 0.59 kPa to 2.45 kPa and that 
the terminal hepatic venule pressure was 0.49 kPa. As we shall demonstrate, this ap- 
plied pressure differential results in a convective flow level that is determined primarily 
by the effective permeability of the lobule. Thus various liver damage scenarios can be 
expected to affect this flow. This aspect of the modelling is of practical importance and 
will be explored in more detail in a separate publication. 
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Table 2 Lobule regular lattice flow parameters 



Parameter 


Characteristic (SI) unit 


STARS unit 


Sinusoid Grid Cell Size 


6 urn 


0.0006 cm 


Sinusoid Porosity: cp s]n 


0.7854 


0.7854 


Sinusoid Permeability: K s]n 


1.125 Mm 2 


1.140 Darcy 


Sinusoid Effective Diffusion: D sjn 


4.2 x 10~ 10 m 2 /sec 


2.5 x 10~ 4 cm 2 /min 


Tissue Grid Cell Size: 2a 


24 urn 


0.0024 cm 


Tissue Porosity (Ideal): cp t[S 


0.4764 


0.4764 


Tissue Permeability (Ideal): K t]s 


1.230 Mm 2 


1.246 Darcy 


Tissue Porosity (Base): cp t]s 


0.2382 


0.2382 


Tissue Permeability (Base): K t]s 


7.35 x 10~ 2 Mm 2 


7.45 x 10" 2 Darcy 


Tissue Porosity (ECM): <p tis 


0.1191 


0.1192 


Tissue Permeability (ECM): K tis 


6.883 X10" 3 Mm 2 


6.97 x 10" 3 Darcy 


Tissue Effective Diffusion: D tis 


4.2 x 10~ 11 m 2 /sec 


2.5 x 10" 5 cm 2 /min 


Blood Viscosity: fj 


3.5 x 10" 3 Pa-sec 


3.5 cpoise 



(1 Darcy = 0.9869 um 2 in engineering permeability units). 



For multicomponent flow, the model tracks the compositions (molar or mass 
fractions) of all components in the fluid. In this paper, we focus on the injected 
drug paclitaxol (PAC) and the Phase I transformed metabolite 6-hydroxypaclitaxel 
(PAC-OH). In addition to convective transport a diffusive flux contribution of 

jf k =D ik V k (p Xi ) (7) 

is considered with ]f k being the molar diffusive flux and D ik the diffusion constant 
of species i in direction k. The estimated diffusion constant in all directions used 
here is based on a molecular weight rescaling of glucose diffusion. Here, glucose 
diffusion coefficient in water is taken as a basic reference value for comparison, 
(E. L. Cussler, "Diffusion: Mass Transfer in Fluid Systems")* 

9 

£^ = 7.1xl0- 10 (8) 

5 sec 

With a cubic root of the molecular weight ratio of glucose to PAC used as conversion 
factor (Factor = (180/854) 033 = (1/4.74) 033 = 1/1.68), an estimated effective diffusion 
constant for PAC is 

7.1 x 1(T 10 in m 2 

Dpact = — = 4.2 x 10- 10 9 

1.68 sec 

Tissue effective diffusion value should be less; here we employ an order of magnitude 
reduction in the value of D 

7.1 x KT 11 n m 2 

Dpact = — = 4.2 x 10 " n 10 

1.68 sec 

Effective diffusion constants for PAC-OH are assumed identical to PAC values. 
These values are converted to the simulation units of cm 2 /min and also summa- 
rized in Table 2. 

The drug paclitaxel was used as a reactive tracer, and its Phase I metabolism 
was modeled using the general formula of one paclitaxel (PAC) molecule being 
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transformed into the metabolite (PAC-OH) by the cytochrome P450 (CYP) 
isozyme CYP2C8 [23]. 

PAC + CYP > PAC-OH + CYP (11) 

The enzyme only exists in grid cells containing hepatocytes so all reaction is localized 
in these sites. In this paper, saturable Michaelis-Menten kinetics are assumed, defined 
by a maximum rate v max (in units of molar fraction/min) and half saturation value K m 
(in units of molar fraction). This reaction proceeds in a linear manner at a rate charac- 
terized by k = v mSLX /K m (in units of min 1 ) when injection concentrations are much 
below the half saturation value. 

Reaction parameter values are based on the work of Vaclavikova et al. [24] who 
measured directly PAC conversion to PAC-OH kinetics without any tissue distribu- 
tion issues and uptake by the cell itself, as they use microsomes directly as the 
source of CYP. As such, any bottlenecks associated with drug uptake should imply 
that reaction rates would be slower than those based on parameters values given by 
Vaclavikova et al. [24]. We are modeling tissue distribution effects separately based 
on our lobule model. Table 3 summarizes the reaction parameters. 

The simulations were performed using the STARS advanced process simulator 
designed by the Computer Modelling Group (CMG) Ltd. in Calgary, Alberta, to model 
the flow and reactions of multiphase, multicomponent fluids through porous media 
[25-27]. Additionally, STARS has earlier been used to model reactive flow processes in 
cortical bone [28-31] as well as through the intervertebral disk [32]. 

Results 

Non-reactive flow characteristics 

As discussed above, flow is induced on the regular lattice of Figure lb by applying 
a pressure difference across the inlet and outlet points. With the chosen lobule flow 
parameters for porosity, permeability, and blood viscosity, this translates to a steady 
flow rate 2.1 cm 3 /min as illustrated in Figure 2a. A short timescale of about 
2.0 x 10~ 5 min needed to establish this pressure gradient is also illustrated in 
Figure 2b by expanding the time axis. Figure 3 demonstrates the steady state 
velocity profile throughout the lattice, illustrating both the diverging/ converging 
nature of the flow near the inlet and outlet points at the top-left and bottom-right 
portion of the grid (ie injector and producer, respectively), as well as the orders of 
magnitude difference of the flows in the sinusoids and tissues, respectively. (This 
plot uses a logarithmic colour scale axis). 

When blood with a relative composition of 1 micro-gram paclitaxel (1.8 x 10 s 
mole fraction) is infused into the lattice assuming nonreactive hepatocytes, the time 



Table 3 Paclitaxel kinetic elimination Michael is-Menten parameters (converted* from 
Vaclavikova et al., their Table 4) 



Parameter 


Characteristic (SI) unit 


STARS unit 


Maximum rate v max 


0.06 uM/min 


1 .08 x 1 0" 9 molefrac/min 


Half saturation constant K m 


10.0 uM 


1.8 x 10" 7 molfrac 


Linear rate v max /K m 


6.0 x 10" 3 min" 1 


6.0 x 10" 3 min" 1 



*Their Table 4 quotes i/ max = 61 picomol/mg_protein/min. Using their microsomal protein concentration of 1 mg/ml, 
these numbers convert to v max = 0.06 uM/min. 
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Figure 2 Injected flow (PAC). (a) Steady state flow across the lobule, (b) Short time flow transients. 



required to traverse the lattice is approximately 1 min as demonstrated in Figure 4. 
This production profile is convective flow dominated as the addition of diffusion 
minimally alters the production profile. 

The evolution of the paclitaxel concentration on the lattice followed a spatially 
homogeneous progression (Figure 5), which shows the increasing levels of injected drug 
after 0.01 min and 0.14 min. By 0.14 min, paclitaxel is being seen at the outlet of the 
lobule. After 0.5 min, paclitaxel completely covers the lattice (not shown). 

If the diffusive flow contribution is removed, however, the paclitaxel profiles on the 
lattice are significantly different. As is also illustrated in Figure 5, again at 0.01 min and 
0.14 min, a distinct two-scale behavior is noted, whereby the sinusoids are first infused 
with the drug, and only at later times do the drug levels in the tissue approach injected 
concentration levels. This behavior reflects the convective levels of flow in the sinusoids 
and tissues noted earlier (Figure 3). A further comparison of Figure 5 cases reveals that 
the sinusoid drug concentration levels in the two cases are similar, however, explaining 



[ctn/min] 

_1.36e+1 
B3,55e+0 

l9.31e-1 

2.44e-1 

- 6.40e-2 
— 1 .68e-2 

- 4.39e-3 
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J 3.02e-4 

7.91e-5 
2.07e-5 



Figure 3 Steady state velocity profile across the lobule. Color bar is in cm/min. 
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Figure 4 Non-reactive PAC drug propagation across the lobule, with and without diffusion effects. 
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Figure 5 Non-reactive PAC profiles across the lobule wit and without diffusion, (a) PAC at 0.01 min, 
with diffusion effect, (b) PAC at 0.14 min, with diffusion effect, (c) PAC at 0.01 min, no diffusion effect, 
(d) PAC at 0.14 min, no diffusion effect. Color bar is in molfrac. 
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the similar drug production characteristics note in Figure 4, as paclitaxel is produced 
directly from the sinusoids. 



Base case reactive flows 

The effects of paclitaxel drug metabolism by hepatocytes are next considered. Here the 
base case reaction parameters of Table 3 are employed, and the same injected paclitaxel 
concentration (1.8 x 10" 8 mole fraction) is considered. With the employed reaction half 
saturation constant value of 1.8 x 10" 7 mole fraction, this injection level implies the 
Michaelis-Menten model reduces to an almost linear reaction scheme. 

Figure 6 illustrates injected drug and produced drug and metabolite production for 
this case. Again it is emphasized that both PAC and PAC-OH have assumed equal 
diffusive flow contributions, as these are components of very similar size. Essentially at 
this reaction rate, all injected paclitaxel is converted to metabolite by the lobule 
hepatocytes. The production profile of PAC-OH here is identical to the production pro- 
file of PAC in the non-reaction case, as shown in Figure 4. Figure 7 shows the PAC and 
PAC-OH profiles across the lobule lattice at 0.01 min, 0.14 min, and 0.50 min, respec- 
tively. The PAC concentrations in the sinusoids and the PAC-OH concentrations in the 
tissue are equivalent to the PAC concentrations in both sinusoids and tissue for the 
non-reacted case (Figure 5). Figure 7 also shows most clearly there is an inlet distance 
over which the reaction conversion time is not fast enough to convert the injected 
paclitaxel. 

Figure 8 illustrates injected drug and produced drug and metabolite production for 
the same case except that diffusive transport has been removed. In contrast to Figure 6 
with diffusion, there is now only a limited amount of conversion of PAC to PAC-OH 
even at long times. The PAC and PAC-OH profiles at various times (0.01 min, 
0.14 min, and 0.50 min) as shown in Figure 9, confirm this behavior where it is shown 
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Figure 6 Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and 
base case metabolism. 
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that the PAC concentration in the sinusoids propagates throughout the lobule, while 
the PAC-OH concentrations in the tissue increase less rapidly and up to a lower level 



Reactive flow sensitivities 

In this section, the consequences of the chosen reaction parameters are illustrated. 
Figure 10 shows production behavior with a 100-fold reduction in maximum reaction 
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Figure 8 Reactive PAC and PAC-OH drug propagation across the lobule, without diffusion effects 
and base case metabolism. 
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rate and with diffusion effects included. The small level of produced PAC indicates that 
reaction rates must be reduced to about this level before any significant change in drug 
production behavior can be expected. Drug distribution in the lobule for this case is 
shown in Figure 11 for the times 0.01 min, 0.14 min, and 0.50 min. This figure should 
be contrasted with Figure 9. Here the early time results and upstream results for PAC 
distributions at longer times are quite different, reflecting the reduced reaction rate. 
However, the later time and downstream results for PAC-OH distribution resemble 
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Figure 10 Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and 
100-fold reduced metabolism. 
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quite closely the faster reaction limit. Here the propagation time across the lobule gives 
enough time to compensate for changes in reaction rate. In summary, reaction rates 
larger than the base case or even 10-fold reduction from base case can be expected to 
produce very similar drug production behavior and differences only in the inlet region 
of the lobule are to be envisioned. 

In contrast, once a critical time-scale is crossed, much more significant changes in 
drug distribution behavior can be expected, both internally throughout the lobule and 
in terms of produced profiles. Figure 12 shows drug metabolite production behavior 
with a 1000-fold reduced metabolic rate, and including diffusive mixing. Here almost 
equal levels of PAC and PAC-OH are seen exiting the lobule. At 0.01 min almost no 
PAC-OH is converted in the lobule at this rate (see Figure 13), while at later times 
(0.14 min and 0.50 min), converted PAC-OH starts to be seen at the outlet regions at 
levels similar to PAC. Essentially, the inlet region behavior occurring at faster reaction 
rates now covers the whole lobule region. 

Finally, sensitivities to injected PAC concentrations were explored, utilizing injection 
concentrations of 1.8 x 10" 7 and 1.8 x 10" 6 mole fractions (i.e. clearly above that of the 
base case 1.8 x 10" 8 ). In these runs, the base case reaction parameters were maintained. 
In particular, the half saturation value of 1.8 x 10" 7 was employed, indicating that the 
linear, intermediate, and saturation levels of the Michaelis-Menten expression were 
being probed with the three injected concentration levels. As illustrated in Figure 14 
for the runs without diffusion, the production profiles of PAC-OH remained unchanged 
for each case, as long as the production maxima were rescaled to the corresponding 
injection concentrations. Apparently, with fast reaction rates, the Michaelis-Menten 
form had little impact on production behavior. 
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Figure 12 Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and 
1000-fold reduced metabolism. 



Analysis 

The results we have presented can be rationalized by a comparison of process time- 
scales. Calculation of breakthrough times can be based on two concepts: either only 
the sinusoids are accessible or the whole lobule (tissue + sinusoid) is accessible to 
injected species. The sinusoid pore volume in our element is 4.77 x 10~ 7 cm 3 while the 
complete lobule element volume is 7.344 x 10" 7 cm 3 . Because the steady state flow rate 
in our model is 2.44 x 10~ 6 cm 3 /min (see Figure 2), this means the breakthrough time is 




Figure 13 Reactive PAC and PAC-OH profiles across the lobule with diffusion effects and 1000-fold 
reduced metabolism, (a) PAC at 0.01 min, (b) PAC at 0.14 min, (c) PAC at 0.50 min, (d) PAC-OH at 
0.01 min, (e) PAC-OH at 0.14 min, (f) PAC-OH at 0.50 min. Color bar is in molefrac. 
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Figure 14 PAC-OH metabolite production levels from various injected PAC concentrations. 



4.77/24.4 = 0.20 min for just sinusoid accessibility and 7.34/24.4 = 0.30 min for the 
whole lobule-sampled space. These should be viewed as two limiting vertical lines on 
the time history plots as two ideal limits without any diffusion or mixing effects 
(physical or numerical). It is noted for example that our time history plot of PAC pro- 
duction with no reaction and with or without physical diffusion (see Figure 4) has a 
produced concentration of 0.9 x 10" 8 (i.e. half of the injected 1.8 x 10" 8 concentration) 
at 0.19 min, about what is expected. The main point here is that most of the pro- 
duction behavior differences for our various cases should lie between these two ideal 
"half-value" limits. 
The next timescale is governed by a "pressure diffusion" coefficient 



K 



pres 



<pf*cf 



(12) 



Here Cf is an effective compressibility accounting for both fluid and tissue structure 
effects. Fluid (water compressibility) is of the order of 5 x 10" 7 kPa" 1 . For liver (soft 
tissue) structural compressibility, we have chosen 1.8 x 10" 5 kPa" 1 . Using these choices 
and the base case parameters of this study from Table 2, we obtain 



D pres =l.5 x 10 4 



cm" 
min 



(13) 



This parameter essentially describes the time taken for pressure to come to a steady 
state distribution as follows. Utilizing a characteristic distance d = 0.15 cm (the lobule 
element size for the pressure calculation), this time is then 



D„ 



= 1.5 x 10" 6 min 



(14) 



Figure 2 also illustrates this characteristic time. This timescale is essentially a function 
of fluid properties and lobule structure (through 0 and K). If we were to consider 
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pulsatile flow effects caused by hepatic artery inflow, this timescale would be much 
more important to the general process description and a more precise definition of 
compressibility might be warranted. For the present, these numbers just indicate that 
steady state pressure is achieved more quickly than other process effects. 

The third timescale is determined by particle diffusion. Here we have chosen 
diffusion constants based on paclitaxel size and simple estimates of tortuosity. The pa- 
rameter choices used here are D sin = 2.5 x 10~ 4 cm 2 /min and D tis = 2.5 x 10~ 5 cm 2 /min. 
As seen from our simulations, for well-perfused tissue convective effects operate 
globally over the whole domain while diffusion smooths concentration profiles locally. 
Thus with a choice of characteristic distance d = 10 \im = 1 x 10~ 3 cm, the times 
required for particles to diffuse are 

d 2 

^j(sin) = 7^=4.0 x 1(T 3 min (15) 

"sin 

d 2 

T diff{ tis) = 7^=4.0 x 10- 2 min (16) 

Our diffusion values should be viewed as highly optimistic. In particular, pact- 
litaxol is normally not molecularly dissolved, but rather it is some type of micellar 
complex with Cremophor EL surfactant, so the effective diffusion coefficient for this 
complex is probably one or more orders of magnitude smaller than what has been 
estimated. Thus the limits of diffusion and non-diffusion cases are meaningful 
extremes of what might be expected, for small molecules and large nanoparticles, 
respectively. 

The final timescale is reaction rate. The base limiting reaction rate 6 x 10 3 min" 1 , also 
utilizing a model tissue volume to bulk volume scale factor of SF ~ 0.76 x 10~ 6 , converts 
to a reaction time (reaction half-life) of 

SF ln(2) 4 , v 

T reac = V=1.0 x l(T 4 min (17) 

6 x 1(T 3 V 7 

This is essentially seen as more rapid than or comparable to the other timescales con- 
sidered (see the early time PAC-OH tissue concentration level appearance in Figures 7d 
and 9d). Reducing this basic reaction rate by a factor of 100 or 1000 causes the reaction 
process to be more similar to the other timescales and different production profiles of 
PAC and PAC-OH result, as has been shown. Table 4 summarizes the relevant assumed 
timescales. 

Table 4 Assumed lobule process time constants 

Process Time 

Convective transit time (sinusoid network only) 0.200 min 

Pressure relaxation time Constant (in sinusoids) 1.5 x 10" 5 min 

Diffusion relaxation time constant (CYP/CYP-OH in sinusoids) 4.0 x 10~ 3 min 

Diffusion relaxation time constant (CYP/CYP-OH in tissue) 4.0 x 10" 2 min 

Base case metabolic uptake/elimination time constant 1 .0 x 1 0" 4 min 
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Conclusion 

In pharmacokinetics, lattice models are introduced to address the non-heterogeneity of 
the organs on the drug distribution that has a significant impact on drug propagation 
throughout the body as shown by analysis of clinical data. 

Here we utilize the interpretation of the liver as an ensemble of islands of metabolic 
activity and focus on the liver lobule itself. In contrast to most of earlier studies that 
assumed that drug molecules only randomly propagate (ie diffuse) through the system, 
we have emphasized the dominant role of convection in well-vascularized tissue with a 
given structure. We have utilized an idealized representation to analyze the factors 
affecting drug propagation and metabolism. The lobule is divided into hepatocyte cells 
that are interlaced with narrower sinusoidal grid cells. These cells are connected by 
constant permeability throughout the entire system. The drug molecules convectively 
flow through the sinusoidal along with blood (water here) and diffuse to the hepatocyte 
where metabolisms are taking place. A sensitivity analysis of convective, diffusive and 
reaction parameters are performed and estimates of their contributions are presented. 
We have compared the drug concentration levels observed exiting the lobule with their 
predicted detailed distribution inside the lobule, assuming that most often the former is 
accessible information while the latter is not. As such, we establish how traditional 
pharmacokinetic analysis might be reflective of the spatial distribution of the drug in 
the lobule, and situations when this might not hold. Interestingly, our network models 
including dispersive effects often correspond to the "well-stirred" compartment models, 
such that relatively uniform steady-state concentration levels occur throughout the 
lobule (if one ignores the smaller inlet mixing zone). Conversely, simulations on our 
network models without explicit dispersive mixing often correspond to modified 
"parallel tube" models, such that observed concentration profiles change along the 
length of the tubes (ie sinusoids). Here our modified tube network structure allows 
cross sinusoids as well. These comments reflect Figures 3, 5, and 9. 

This is the first paper of series of papers on physiologically-based lattice models for 
liver. In this paper, we consider an idealized lobule lattice in order to understand the 
basic functionality of the unit and underlying mechanisms through simulations and also 
to set a basis for future studies. In following papers, we will expand the analysis to 
include drug propagation & metabolic sensitivities associated with variations in lobule 
structure, which could reflect extents of liver damage, and how our modelling approach 
might be used to generate flows on realistic, reconstructed images of lobule structure. 
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